# =================================================================================
# remove everything
rm(list=ls())

# ---------------------------------------------------------------------
# Drop outliers in dataset "Data" according to the "v"th column
# at each percentage "p" at bottom and top. The default value of
# p is 0.005. v is a vector of numbers.
# ---------------------------------------------------------------------

outlier <- function(Data, v, p=0.005) {
  out <- c()
  for(i in v) {
    out <- rbind(out, quantile(Data[,i], probs=c(p, 1-p), na.rm = TRUE))
  }
  j <- 1
  for(i in v) {
    Data <- Data[Data[,i]>out[j,1] & Data[,i]<out[j,2], ]; print(dim(Data))
    j <- j+1
  }
  return(Data)
}
# ---------------------------------------------------------------------

# Load packages
library(foreign)
library(plm) # ---------------

# Load data
Data <- read.dta('Ind39.dta')

Data <- outlier(Data, which(names(Data) %in% c("sm", "Employment", "TotalProfit", "ROA", "Lev", "TotalAssets")), p=0.005)
Data[,which(names(Data)=="Sub")] <- (Data$Sub>0)

# Delete obs according to variable "Foreign"
Data <- subset(Data, F1>=0 & F1<=1 & !is.na(Province) & sm<=1 & !is.na(Employment) & !is.na(age) & SOE >=0 & SOE <=1 & ExpR >=0 & ExpR <=1 & age >0 & age <500 & TotalAssets>0)

# -------------------------------
# # Generate Propensity Score
ConV <- cbind(log(Data$Employment), log(Data$age), Data$SOE, (Data$Subsidy), Data$ExpR, scale(Data$TotalProfit), Data$ROA, Data$Lev, log(Data$TotalAssets), Data$coast, Data$T)

N.col <- ncol(ConV); 
F.beta <- function(a) {
  XX <- cbind(rep(1, nrow(ConV)), ConV)
  a1 <- 1; a2 <- (sum(Data$F1==0)/nrow(ConV)*100)/(exp(ConV%*%a)+1)
  obj <- -sum(log(dbeta(Data$F1, a1, a2)+1))/nrow(ConV)
  return(obj)
}

set.seed(123)
initial <- -abs(rnorm(N.col, 0, 0.1))
mod.p <- optim(initial, F.beta, method = "BFGS")

PS.d <- dbeta(Data$F1, 1, (sum(Data$F1==0)/nrow(ConV)*100)/(exp(ConV%*%mod.p$par)+1) )

F.beta2 <- function(a) {
  a1 <- 1; a2 <- (sum(Data$F1==0)/nrow(ConV)*100)/(exp(a)+1)
  obj <- -sum(log(dbeta(Data$F1, a1, a2)+1))/nrow(ConV)
  return(obj)
}
mod.p2 <- optim(3, F.beta2, method = "BFGS")

PS.num <- dbeta(Data$F1, 1, (sum(Data$F1==0)/nrow(ConV)*100)/(exp(mod.p2$par)+1) )
PS <- as.vector(PS.num/PS.d)
Data$PS <- PS

# Generate weights & variables related to weights
for(i in 1:nrow(Data)) {
  s1.t <- (Data$Year==Data$Year[i])*(Data$Province==Data$Province[i]); s1.t[i] <- 0
  s1.b <- sum(s1.t)
  s1 <- if(s1.b>0) s1.t/s1.b else rep(0, length(s1.t));
  Data$wk1[i] <- sum(s1*Data$k1); Data$wm1[i] <- sum(s1*Data$m1, na.rm = TRUE); Data$wl1[i] <- sum(s1*Data$l1, na.rm = TRUE)
}

# Generate year dummies
Data$TD <- model.matrix(~as.factor(Data$Year)-1); Data$TD1 <- model.matrix(~as.factor(Data$Year-1)-1)[,-1]

# Calculate the length of loop
firm.list <- unique(Data$Firm)
R <- ceiling(length(firm.list)/50)

# empty obj.
JKpw <- JKpF <- JKpsw <- JKptil <- JKbeta <- JKcoef <- c()

for(ii in 1:R) {
  
  print(ii)
  # get jackknife samples
  from <- 50*(ii-1)+1; to <- 50*ii
  index <- Data$Firm %in% firm.list[from:to]
  D.jk <- Data[!index, ]
  
  # =============================================
  # Step One
  h.lnbe <- mean(log(D.jk$sm)); h.eta <- - log(D.jk$sm) + h.lnbe; h.e <- mean(exp(h.eta))
  h.betam <- exp(h.lnbe)/h.e

  # =============================================
  # Step Two
  # Define functions: "wphi" and "phi"
  wphi <- function(beta) {
    (1-h.betam)*D.jk$wm1 - h.lnbe - log(D.jk$ppi1/D.jk$ppii1) - beta[1]*D.jk$wk1 - beta[2]*D.jk$wl1 - D.jk$TD1%*%beta[3:10] 
  }
  
  phi <- function(beta) {
    (1-h.betam)*D.jk$m1 - h.lnbe - log(D.jk$ppi1/D.jk$ppii1) - beta[1]*D.jk$k1 - beta[2]*D.jk$l1 - D.jk$TD1%*%beta[3:10]
  }

  # -----------------------------------------------
  # Optimization
  objfn <- function(beta){
    y.star <- D.jk$y - h.betam*D.jk$m - beta[1]*D.jk$k - beta[2]*D.jk$l - D.jk$TD%*%beta[3:11]
    p1 <- as.vector(phi(beta)); p2 <- as.vector(wphi(beta)); 
    var.poly <- poly(cbind(p1, D.jk$F1, p2), degree = 2, raw = TRUE)
    mod1 <- lm(as.vector(y.star)~var.poly, weights = D.jk$PS)
    OBJ <- sum(mod1$residuals^2)
    return(OBJ)
  }
  
  set.seed(1)
  initial <- c(0.05, 0.1, 0.05, 0,0,0,0,0,0,0,0,0)
  target <- optim(initial, objfn, method = "Nelder-Mead")
  
  # results - elas
  h.betak <- target$par[1]; h.betal <- target$par[2]
  
  G.coef <- function(beta){
    y.star <- D.jk$y - h.betam*D.jk$m - beta[1]*D.jk$k - beta[2]*D.jk$l - D.jk$TD%*%beta[3:11]
    p1 <- as.vector(phi(beta)); p2 <- as.vector(wphi(beta)); 
    var.poly <- poly(cbind(p1, D.jk$F1, p2), degree = 2, raw = TRUE)
    mod1 <- lm(as.vector(y.star)~var.poly, weights = D.jk$PS)
    return(mod1$coef)
  }
  # results - partial effects
  coef <- G.coef(target$par)[2:10]
  
  # Define functions: "wphi" and "phi" for the full sample
  wphi <- function(beta) {
    (1-h.betam)*Data$wm1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$wk1 - beta[2]*Data$wl1 - Data$TD1%*%beta[3:10] 
  }
  
  phi <- function(beta) {
    (1-h.betam)*Data$m1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$k1 - beta[2]*Data$l1 - Data$TD1%*%beta[3:10]
  }
  
  w1 <- as.vector(phi(target$par)); sw1 <- as.vector(wphi(target$par))
  
  p.w  <- coef[1] + 2*coef[2]*w1 + coef[4]*Data$F1 + coef[7]*sw1
  p.F  <- coef[3] + coef[4]*w1 + 2*coef[5]*Data$F1 + coef[8]*sw1
  p.sw <- coef[6] + coef[7]*w1 + coef[8]*Data$F1 + 2*coef[9]*sw1
  
  # TIL
  D <- data.frame(Data[,c('Firm','Year','Province')], pw=p.w, pF=p.F, psw=p.sw)
  D <- pdata.frame(D, index = c("Firm", "Year"))
  D$pF1 <- lag(D$pF, 1); 
  
  D$temp1 <- ave(D$pF1, D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
  D$temp1[!is.na(D$pF1)] <- D$temp1[!is.na(D$pF1)] - D$pF1[!is.na(D$pF1)]
  D$temp2 <- ave(D$pF1, D$Year, D$Province, FUN = function(x) sum(!is.na(x)))
  D$temp2[!is.na(D$pF1)] <- D$temp2[!is.na(D$pF1)] - 1
  D$sp=D$temp1/D$temp2
  p.til <- as.numeric(D$psw*D$sp); summary(p.til)
  
  # Store bootstrap results
  JKpw <- cbind(JKpw, p.w); JKpF <- cbind(JKpF, p.F); JKpsw <- cbind(JKpsw, p.sw); JKptil <- cbind(JKptil, p.til)
  JKbeta <- cbind(JKbeta, c(h.betak, h.betal, h.betam))
  JKcoef <- cbind(JKcoef, coef)
  
}

apply(JKbeta, MARGIN = 1, sd)
summary(JKpw[,1]); apply(apply(JKpw, MARGIN = 2, summary), MARGIN = 1, sd)
summary(JKpF[,1]); apply(apply(JKpF, MARGIN = 2, summary), MARGIN = 1, sd)
summary(JKpsw[,1]); apply(apply(JKpsw, MARGIN = 2, summary), MARGIN = 1, sd)
summary(JKptil[,1]); apply(apply(JKptil, MARGIN = 2, summary), MARGIN = 1, sd)

save(JKpF, JKpw, JKpsw, JKptil, JKbeta, JKcoef, file = "Jack_s_Dt_IPW.RData")
